Using FISSA with Suite2p

suite2p is a blind source separation toolbox for cell detection and signal extraction.

Here we illustrate how to use suite2p to detect cell locations, and then use FISSA to remove neuropil signals from the ROI signals.

The suite2p parts of this tutorial are based on their Jupyter notebook example.

Note that the below results are not representative of either suite2p or FISSA performance, as we are using a very small example dataset.

Reference: Pachitariu, M., Stringer, C., Dipoppa, M., Schröder, S., Rossi, L. F., Dalgleish, H., Carandini, M. & Harris, K. D. (2017). Suite2p: beyond 10,000 neurons with standard two-photon microscopy. bioRxiv: 061507; doi: 10.1101/061507.

Imports

! Important note: if you want to run Suite2p and FISSA in the same python instance (as we do in this notebook), you have to make sure multiprocessing pools are started using the 'spawn' context, by having the following on top of your code: See here for more information: https://docs.python.org/3/library/multiprocessing.html

from multiprocessing import set_start_method
set_start_method("spawn") # this is necessary to make FISSA's multiprocessing play nicely with Suite2p's

# FISSA toolbox
import fissa

# suite2p toolbox
import suite2p

# numpy toolbox
import numpy as np

# Plotting toolbox, with notebook embedding options
import holoviews as hv
%load_ext holoviews.ipython
%output widgets='embed'

Run suite2p

# Set your options for running
ops = suite2p.default_ops()  # populates ops with the default options

# Provide an h5 path in 'h5py' or a tiff path in 'data_path'
# db overwrites any ops (allows for experiment specific settings)
db = {
    'h5py': [],  # a single h5 file path
    'h5py_key': 'data',
    'look_one_level_down': False,  # whether to look in ALL subfolders when searching for tiffs
    'data_path': ['exampleData/20150529'],  # a list of folders with tiffs 
                                            # (or folder of folders with tiffs if look_one_level_down is True,
                                            # or subfolders is not empty)
    'save_path0': './',  # save path                                      
    'subfolders': [],  # choose subfolders of 'data_path' to look in (optional)
    'fast_disk': './',  # string which specifies where the binary file will be stored (should be an SSD)
    'reg_tif': True,  # save the motion corrected tiffs
    'tau': 0.7,  # timescale of gcamp6f
    'fs': 1,  # sampling rate
    'spatial_scale': 140, # scale of cells
    'batch_size': 288 # size of each tif
}

# Run one experiment
opsEnd = suite2p.run_s2p(ops=ops, db=db)

{'h5py': [], 'h5py_key': 'data', 'look_one_level_down': False, 'data_path': ['exampleData/20150529'], 'save_path0': './', 'subfolders': [], 'fast_disk': './', 'reg_tif': True, 'tau': 0.7, 'fs': 1, 'spatial_scale': 140, 'batch_size': 288}
FOUND BINARIES AND OPS IN ['./suite2p/plane0/ops.npy']
>>>>>>>>>>>>>>>>>>>>> PLANE 0 <<<<<<<<<<<<<<<<<<<<<<
NOTE: not running registration, plane already registered
binary path: ./suite2p/plane0/data.bin
NOTE: Applying builtin classifier at /home/sander/anaconda3/lib/python3.8/site-packages/suite2p/classifiers/classifier.npy
----------- ROI DETECTION
Binning movie in chunks of length 01
Binned movie [500,172,152] in 0.13 sec.
NOTE: FORCED spatial scale ~48 pixels, time epochs 1.00, threshold 20.00 
0 ROIs, score=237.58
Detected 23 ROIs, 0.68 sec
After removing overlaps, 23 ROIs remain
----------- Total 0.87 sec.
----------- EXTRACTION
Masks created, 0.06 sec.

/home/sander/anaconda3/lib/python3.8/site-packages/suite2p/extraction/extract.py:111: NumbaTypeSafetyWarning: unsafe cast from uint64 to int64. Precision may be lost.
  Fi[n] = np.dot(data[:, cell_ipix[n]], cell_lam[n])

Extracted fluorescence from 23 ROIs in 864 frames, 3.09 sec.
----------- Total 3.16 sec.
----------- CLASSIFICATION
['npix_norm', 'skew', 'compact']
----------- Total 0.03 sec.
----------- SPIKE DECONVOLUTION
----------- Total 0.00 sec.
Plane 0 processed in 4.06 sec (can open in GUI).
total = 4.07 sec.
TOTAL RUNTIME 4.07 sec

Load the relevant data from the analysis

# Extract the motion corrected tiffs (make sure that the reg_tif option is set to true, see above)
images = './suite2p/plane0/reg_tif'

# Load the detected regions of interest
stat = np.load('./suite2p/plane0/stat.npy', allow_pickle=True)  # cell stats
ops = np.load('./suite2p/plane0/ops.npy', allow_pickle=True).item()
iscell = np.load('./suite2p/plane0/iscell.npy', allow_pickle=True)[:, 0]

# Get image size
Lx = ops['Lx']
Ly = ops['Ly']

# Get the cell ids
ncells = len(stat)
cell_ids = np.arange(ncells)  # assign each cell an ID, starting from 0.
cell_ids = cell_ids[iscell==1]  # only take the ROIs that are actually cells.
num_rois = len(cell_ids)

# Generate ROI masks in a format usable by FISSA (in this case, a list of masks)
rois = [np.zeros((Ly, Lx), dtype=bool) for n in range(num_rois)]

for i, n in enumerate(cell_ids):
    # i is the position in cell_ids, and n is the actual cell number
    ypix = stat[n]['ypix'][~stat[n]['overlap']]
    xpix = stat[n]['xpix'][~stat[n]['overlap']]
    rois[i][ypix, xpix] = 1

Run FISSA with the defined ROIs and data

output_folder = 'fissa_suite2p_example'
exp = fissa.Experiment(images, [rois[:ncells]], output_folder)
exp.separate()

Reloading previously prepared data...
Reloading previously separated data...
Doing region growing and data extraction....
Doing signal separation....

/home/sander/anaconda3/lib/python3.8/site-packages/numpy/core/_asarray.py:136: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order, subok=True)

Plot the resulting ROI signals

%%opts Curve {+axiswise}


def plot_cell_regions(roi_polys, plot_neuropil=False):
    '''
    Plot a single cell region, using holoviews.
    '''
    out = hv.Overlay()

    if plot_neuropil:
        # Plot the neuropil as well as the ROI
        n_region = len(roi_polys)
    else:
        # Just plot the ROI, not the neuropil
        n_region = 1

    for i_region in range(n_region):
        for part in roi_polys[i_region]:
            x = part[:, 1]
            y = part[:, 0]
            out *= hv.Curve(zip(x, y)).opts(color='w')

    return out


i_trial = 0

# Generate plots for all detected regions
region_plots = {
    i_cell: plot_cell_regions(exp.roi_polys[i_cell][i_trial])
    for i_cell in range(exp.nCell)
}

# Generate plots for raw extracts and neuropil removed
traces_plots = {
    i_cell: hv.Curve(exp.raw[i_cell][i_trial][0, :], label='suite2p') *
            hv.Curve(exp.result[i_cell][i_trial][0, :], label='FISSA')
    for i_cell in range(exp.nCell)
}

# Generate average image
avg_img = hv.Raster(exp.means[i_trial])

# Generate overlay plot showing each cell location
cell_locs = hv.Overlay()
for c in range(exp.nCell):
    roi_poly = exp.roi_polys[c][i_trial][0][0]
    x = roi_poly[:, 1]
    y = roi_poly[:, 0]
    cell_locs *= hv.Curve(zip(x, y))

# Render holoviews
avg_img * cell_locs * hv.HoloMap(region_plots, kdims=['Cell']) + hv.HoloMap(traces_plots, kdims=['Cell'])

Note that with the above settings for suite2p it seems to have detected more small local axon signals, instead of cells. This can possibly be improved with manual curation and suite2p setting changes, but as noted above these results should not be seen as indicative for either suite2p or FISSA due to the small dataset size.

Also note that the above Suite2P traces are done without suite2p's own neuropil removal algorithm.